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Abstract. Two methods to approximate infinitely divisible random fields are pre- 
sented. The methods are based on approximating the kernel function in the spectral 
representation of such fields, leading to numerical integration of the respective inte- 
grals. Error bounds for the approximation error are derived and the approximations 
are used to simulate certain classes of infinitely divisible random fields. 



1. Introduction 

In many cases, the normal distribution is a reasonable model for real phenomena. 
If one considers the cumulative outcome of a great amount of influence factors, the 
normal distribution assumption can be justified by the Central Limit Theorem which 
states that the sum of a large number of independent and identically distributed 
random variables can be approximated by a normal distribution if the variance of 
these variables is finite. However, many real phenomena exhibit rather heavy tails. 
Stable distributions remedy this drawback by still being the limit distribution of a 
sum of independent and identically distributed random variables, but allowing for an 
infinite variance and heavy tails. 

Stable distributions are a prominent example of the class of infinitely divisible dis- 
tributions which we particularly concentrate on in this paper. Infinitely divisible 
distributions are distributions whose probability measure P is equal to the n-fold 
convolution of a probability measure P n for any positive integer n. The class of 
infinitely divisible distributions comprises further well-known examples such as the 
Poisson, geometric, negative binomial, exponential, and gamma distribution, see |15| . 
p. (page) 21. These distributions are widely used in practice, for instance in finance 
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to model the returns of stocks or in insurance to model the claim amounts and the 
number of claims of an insurance portfolio. 

In order to include time dependencies or the spatial structure of real phenomena, 
random processes may be an appropriate model. If the dimension of the index set 
of the random process is greater than one, random processes are also called random 
fields. An example of infinitely divisible processes are Levy processes which have been 
extensively studied in the literature. 

In this paper, we consider random fields that can be represented as a stochastic inte- 
gral of a deterministic kernel as integrand and an infinitely divisible random measure 
as integrator. The kernel basically determines the dependence structure, whereas the 
infinitely divisible random measure inhibits the probabilistic characteristics of the 
random field. As already noted by pj, practitioners have to try a variety of kernels 
and infinitely divisible random measures to find the model that best fits their needs. 

Once the model is fixed, it is desirable to be able to perform simulations of the 
considered random field. There are several papers that are devoted to this problem. 
In [I], [17] and [22], the fast Fourier transform is used for the simulation of linear 
fractional stable processes, whereas in [6], a wavelet representation of a certain type 
of fractional stable processes was applied to simulate sample paths. Furthermore, [3] 
gives a general framework for the simulation of fractional fields. 

In this paper, we consider infinitely divisible random fields for which the kernel func- 
tions are assumed to be Holder-continuous or bounded which is a less restrictive 
assumption. Based on the respective assumption, we derive estimates for the approx- 
imation error when the kernel functions are approximated by step functions or by 
certain truncated wavelet series. The approximation allows for simulation since the 
integral representation of the random field reduces to a finite sum of random variables 
in this case. 

In Section 02, we present the main results for the approximation error which is made 
when the kernel function is replaced by a step function or a truncated wavelet series. 
Section 0] is devoted to a brief simulation study where we apply the derived formulas 
for the approximation error to the simulation of two particular stable random fields. 
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Finally, in Section [5] we comment on the simulation results and the methods discussed 
in Section 03 



2. Infinitely divisible random fields admitting an integral 

representation 

Let A be an infinitely divisible random measure with control measure A, cf. (com- 
pare) [TB], pp. (page and the following) 455 and [9], pp. 75. Let f t : R d — > R, d > 1, be 
A-integrable for all t £ R 9 , q > 1, that is there exists a sequence of simple functions 
{ft n) }nm, ft n) ■ R d -> R for all t G R 9 , such that 

(a) fV-ft A-a.e., 

(b) for every Borel set B C R d , the sequence {/ ft\x)A(dx)} ne ® converges in 

B 

probability. 



For each t £ M 9 , we define 



/t(z)A(dx) :=plim / / t (n) (x)A(rfx), 

n— >oo J 



cf. [32], p. 460, where plim means convergence in probability, and consider random 
fields of the form 



X(t) = J f t (x)A(dx), t G W. (1) 



Remark 2.1. In |13J. it is shown that X{t) is infinitely divisible for all t G R 9 , cf. The- 
orem 2.7, p. 461 and the Levy form of the characteristic function of an inifinitely 
divisible random variable, p. 456. More generally, it can be shown that the random 
vector (X(ti), X(t n )) is infinitely divisible for all t\, t n G W and n£N, see [TT] . 
Therefore, any random field of the form (pp) is an infinitely divisible random field. 

Example 2.2. Let < a < 2, M be an (independently scattered) a-stable random 
measure on R d with control measure m and skewness intensity ft, see [H], pp. 118. 
Furthermore, we assume that f t G Z/*(R d ,m) if a ^ 1 and f t £ {f £ L 1 (R d , m) : 
J Rd |/(x) In \f(x)\\m(dx) < oo} if a = 1 for all t £ R q . We denote the set of all 
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functions f t satisfying these conditions by T . Then 

X(t) = J f t (x)M(dx), t E R q , 

is an a-stable random field. 

Example 2.3. Let $ be a Poisson random measure on M d with intensity measure 0, 
see pEE], p. 42. Furthermore, we assume that f t is a measurable function on R d for 
each t G M. q . Then we can consider the shot noise field 



X(t)= / ft(x)$(dx), teW, 

which can be written as 



X(t) = Y,M*), (2) 



where \& is the support set of the Poisson random measure, cf. [IB], p. 101. 

Example 2.4. Let Q be a Poisson random measure on (0, oo) x K\ {0} with intensity 
measure fi x v. Here, \i is the Lebesgue measure and v is the Levy measure, cf. [12] and 
[7]. Let G be a Gaussian (2-stable) random measure with Lebesgue control measure 
and skewness intensity {3 = 0. Then 

X(t) = J xl{0 <s< t}Q(ds, dx)-t J xv{dx) + jt + J I{0 < x < t}G(dx) 

R 2 \x\<l R 

is the Levy process with Levy measure u, Gaussian part G and drift 7, where 



7 = E [X 1 - / xv(dx)) 

V J\x\>l J 



\x\>l 

We now consider the cumulant function C*a(a)(^) = ln(Ee l<A ^) of A(A) for a set A in 
the 5-ring A of bounded Borel subsets of M. d which is given by the Levy-Khintchine 
representation 

C M a)(v) = iva(A) - \v 2 b(A) + / (e ivr - 1 - ivrl hhl] (r))U(dr, A), 

* Jr 

where a is a cx-additive set function on A, b is a measure on the Borel a-algebra £>(IR d ), 
and U (dr, A) is a measure on B{R d ) for fixed dr and a Levy measure on B{R) for 
each fixed A e B{R d ), that is U({0},A) = and J R min{l, r 2 }U (dr, A) < 00, cf. [8], 
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p. 605. The measure U is referred to as the generalized Levy measure and (a, b, U) is 
called characteristic triplet. The control measure A can be written as 



X(A) = \a\(A) + b(A) + / min{l, r 2 }U(dr, A), A e A, 



(3) 



where \a\ = a + + a , see [13], p. 456. Furthermore, a and b are absolutely continuous 
with respect to A and we have the formulas 

a(drj) — a(rj)\{drj) , b(drj) = b{rj) \{drj) , U(dr,drj) = V(dr,r})\(drj), 

where V(dr,rj) is a Levy measure for fixed 77, cf. [13], p. 457. 

We now introduce a so-called spot variable L'{rj) with cumulant function 

C L '( v )(v) = iva{rj) - \v 2 b(r,) + / (e ivr - 1 - ivrl { _ ltl] (r))V(dr, 77) 

with 

E(L'(r))) = 0(77)+ / rV(dr,r]), (4) 
Var(L'(r7)) = 6(77)+ / rV(dr, 77) (5) 



if E(L'0?)) and Var(L'(r/)) exist. In [8], p. 607, it is shown that 

c X (t)(v)= [ CvwivfttoMdn). 



We can use the cumulant function of X(t) to obtain the second moment of X(t) (in 
the case it exists): 

E(X(t) 2 )=[ f 2 {y)V^{L'{y))\{dy) + ( [ f t (y)E(L'(y))X(dy)) 2 . (6) 

JR d \JR d J 

As noted by [8], for modelling purposes it is no resriction if we only consider charac- 
teristic triplets (a, b, U) of the form 

a(drj) = a v (rj)u(drj), b(drj) = b v (rj)u(drj), U(dr,drj) = V u (dr,drj)u(drj), 

where v is a nonnegative measure on B(M. d ), d v : M. d — > K and 6^ : M d — > [0, 00) are 
measurable functions, and V u (dr,rj) is a Levy measure for fixed 77. 
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Example 2.5. We choose the Lebesgue measure for v and consider the characteristic 
triplet (a, 0, U) with 

U(dr,drj) = V(dr,r})drj = JL(o t00 )(r)-e~ er drdr) 

a{drj) = a(r))dr) — - (l — e" e ) dr] 
with 9 G (0, oo). Then, by using (J4|) and we get 

mm = \, Var(L'fo)) = 1 
and by (J3j) , the control measure is proportional to the Lebesgue measure with 
^f J ^ fl + 0-28e~ e -e- e fl^V 

= ^ 02 + I ~ e dr ) d ^ 

3. Approximation of infinitely divisible random fields 
We now restrict our setting to the observation window [— T, T] q with T > such that 



X(t) = J f t (x)A(dx), t G [-T, T] q . 



We denote by supp(f t ) the support of f t for each t G [—T,T] q and assume that 

|J supp(f t )c[-A,A] d 
te[-T,T]i 

for an A > 0. Then X(-) can be written as 



X(t)= J f t (x)A(dx), te[-T,T] q . 



[-A,A] d 

Our goal is to approximate sample paths of X for a variety of kernel functions f t , 
t G [—T,T] q . The idea is to approximate the kernel functions f t appropriately such 
that the approximations are of the form 

m(n) 

i=i 

where m(n) G N, <ij G K and g^j : M d — > M is A-integrable. 
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Due to the linearity of the stochastic integral, we get as an approximation X^' of X 

/m(n) 
jT\x)A(dx) = V a t \ g t>i (x)Hdx), t E [— T, T] q . 
L -A,A] d ~[ J[-A,A] d 

If the g t A are simple functions such that 

/ g^(x)A(dx) = Y,9t,i(xi)M&j), i = l,...,m(n), t e [-T,T\«, 

J[-A,A] d j=1 

for some Xj E [— A, A] d , I E N and a partition {Aj} l j=1 of [— A, A] d , then 

m(n) i 

^ (n) W = EE a ^(^)A(A,) 

i=l j=l 

which can be simulated if A(Aj-), j = 1, can be simulated. 

Example 3.1. Let A = M be an a-stable random measure with Lebesgue control 
measure and constant skewness intensity (3. Then 

M(A,)~S a (|A/ /Q ,A0), J = W, 

cf. [H], p. 119, where |Aj| is the volume of Aj and S a (a,f3,0) denotes the stable 
distribution with stable index a, scale parameter a, skewness parameter (3 and lo- 
cation parameter 0. Furthermore, M(Aj), j = 1, are independent since M is 
an independently scattered random measure. A method to simulate a-stable random 
variables can be found in [2]. 

Example 3.2. Let A = $ be a Poisson random measure with intensity measure 6. 
Then 

A(Aj) ~ Poi(Q(Aj)), 

where Poi (Q(Aj)) denotes the Poisson distribution with mean Q(Aj). We note that 
simulating sample paths of X by X^ n ' is not efficient since one can directly exploit 
the structure of X and use 

cf. equation ((2j) in Example 12.31 p. [H 
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Example 3.3. Let A x = Q be a Poisson random measure on (0, oo) xl \ {0} with 
intensity measure fix u, where [i is the Lebesgue measure and v is the Levy measure. 
Let A 2 = G be a Gaussian random measure with Lebesgue control measure and 
skewness intensity (3 = 0. We first approximate the Levy process 

X(t) = J xl{0 <s< t}Q(ds, dx)-t J xv(dx) + -ft + J I{0 < x < t}G(dx) 

R 2 |x|<l R 

by 

K t t 

x K(t) = // ie( ^ I )- 1 / I *, + 7l + / GW 

—K |x|<l 

K t 




xQ(ds,dx) -t J xv(dx) +jt + G([0,t)) 

-K |x|<l 

for some K > 0. We approximate f(x) = x by using a linear combination of some 
simple functions g i7 i = 1, ...,m(n), m(n) G N for all n G N and get for a partition 
{A y }5 =1) I G N, of x [0,t] 

m(n) / 

= Z Z <>i9ifa)Q(&j) - t / + 7* + G([0, t}), 

4=1 ' =1 |X|<1 
for some a%, a m ( n ) G R, where 

g(Aj-) ~Poi(0iXi/)(A i )) and G([0,i\) ~ Af(y/i,0). 

Example 3.4. We choose again the Lebesgue measure for v and the characteristic 
triplet (a, 0, U) from Example 12.51 p. Then 



A(A,)~r(|A 



cf. [H], p. 608, where \ Aj\ is the Lebesgue measure of Aj and r(|Aj|,0) is the gamma 
distribution with probability density function 

/(*) = T^^'"^ W*)- 

In the last formula, T(-) denotes the gamma function. Again, A(Aj), j = 1, ...,n, are 
independent. Due to its distributional property, A is called gamma Levy basis. 



SIMULATION OF INFINITELY DIVISIBLE RANDOM FIELDS 9 

3.1. Measuring the approximation error. 

Approximating the random field X with X^ by taking an approximation ft of the 
kernel functions f t implicitely includes the assumption that X'"' is close to X when 
ft^ is close to ft- We use 

Err s (X(t),X^(t)):= 

to measure the approximation quality of X' n ' for an s > 0. In order to have existance 
of the above integral, we assume from now on that ft, ft e L s ([—A,A] d ,\) for 
all t G [—T,T] q . The goal is then to find a set of functions {/j^jtgRd such that 
Err s (X(t),X ( - n '(t)) is less than a predetermined critical value. 

We see that the problem of approximating the random field X reduces to an approx- 
imation problem of the corresponding kernel functions. 

Let us now consider two special cases where A is an a-stable random measure and a 
Poisson random measure, respectively, to analyse the choice of the error measure. 



x] 



■A,A] C 



\ft 



X 



-fi n) (x)\ s X(dx) 



i 



3.1.1. a-stable random measures. 

Assume that < a < 2 and let M be an a-stable random measure with control 
measure m and skewness intensity (3. Furthermore, if a = 1, assume additionally 
that /3(t) = for all t e [—T,T] q . Consider a set of functions {ft }teR d i where 
ft G F (cf. Example 12.21 p. [3]) for all t G [—T,T} q and n G N. The corresponding 
ct-stable random field is denoted by 

J ft\x)M(dx), t G [-T, T] q . 

[—A.A] d 

We know that IW(i) converges to X{t) in probability if and only if 
j R d \ft(x) — f\: n \x)\ a m(dx) converges to as n goes to infinity, see [HJ, p. 126. 
Therefore, we can use X^ n \t) as an approximation for X(t) if f^ approximates 
ft sufficiently well and X^ n \t) converges to X(t) in probability if and only if 

Err a (X(t),X {n) (t)) -> 0, n -> oo. 

The choice of Err a (X(t), X^^t)) can be further justified as follows. 
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Since X(t) and X^ n \t) are jointly a-stable random variables for all t G [— T, T] 9 , the 
difference X{t)-X^ n \t) is also an a-stable random variable. The scale parameter of 
X(t)-XW(t) is given by 

^( { )4W(<) = ( / \ft( x ) ~ fi n) {x)\ a m{dx)\ , 
cf. [33], p. 122, so that 

Err a (X(t) J X^(t))=a x{t) _^ { ty 

Furthermore, let us consider the quantity 

E\X(t) -X {n) (t)\ p , 0<p<a, 
that is the mean error between X(t) and X^(t) in the D> -sense. 
Since X(t) — X^ n \t) is an a-stable random variable, we have 

E\X(t) - X {n) (t) \ p < oo, < p < a 

and 

E\X(t) -X^ n \t)\ p = oo, p>a. 
For 0<p<a,0<ct<2 and a ^ 1, this quantity can be written as 

E|X(f) - I (n) (t)r) 1/P = c aA {p) ■ a m _ tin){t) , (7) 

and 



where 



2^r(l-§) / ^ a a7T\P/*» /p / «7TN 

l c a,/3 t (p)) = — Foo i . 2 — 3- 1 + P t tan — cos - arctan fj t tan — 

p J r^'sin 11 aw V 2 / \a V 2 / 



L AA]d \ft(x)-fV(x)\ a dx 



i,A] c 

In the case a — 1, equation (jTj) holds if $ = 0, see [H], p. 18. 

We remind that /?(•) is the skewness intensity of the a-stable random measure M. 
The above implies that 

Err a (X(t),X^(t)) = — ^ • (e|X(*) - X™(t)\ p ) 1/P , < p < a. 
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Now assume that a = 1 and f3{t) ^ for at least one t e [— T, T] q . In this case, we 
need to impose an additional condition on the kernel functions in order to guarantee 
the convergence of X^(t) to X(t) in probability. Namely, X^ n \t) converges to X(t) 
in probability if and only if 

S[-A,A] d 1-^^) ~ ft^ ( X )\ m ( dx ) °» 71 -»• 00 

and 

Ji-^-U^) - ft n \x)) In 1/tO*) - ^(sMxMdz) - 0, n 
cf. [H], p. 126. We use 

\ 2 - 

\f t (x)-fT\x)f 2 m(dx) 



oo, 



Err 3/2 (X(t),XW(t)) := 
to measure the approximation error. 



A,A] C 



By using the fact that |xlnx| < m&x{y/x,Xy/x} for x > and that —1 < (3(t) < 1 
for all t G [— T, T] 9 , we have 



< 



< 



< 



AA] d 
A,A] d 

A,A] d 
A,A] C 



(ft(x) - ~fi n \x)) In \f t (x) - ~ft\x)\(3{x)m{dx) 
f t (x) - fT\x))\n\f t (x) - fT\x)\P(x)\m(dx) 
max(|/ t (x) - fi n \x)\ l '\\f t {x) - fi n \x)f 2 )m{dx) 



AA] d 



\f t (x) - ft\x)\ 3/2 m(dx) 



Now, if Err 3 / 2 (X(t),X^ n \t)) tends to as n goes to infinity, then by Lyapunov's 
inequality (see [TB]). 



/ \f t (x) - fi n \x)\*' 2 m(dx) ^ 

J[-A,A]d 

[ \f t {x)-fi n \x)\ l/2 m{dx)^^ n 



n — > oo, 



oo, 



such that 



(ft(x) - ft\x))\n\ft{x) - fr(x)\/3(x)m(dx) - 0, n oo. 



A,A] C 
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Furthermore, once again by using Lyapunov's inequality, 

\ft(x) - fl n) (x)\m(dx) -> 0, n->oo, 

■A,A\ d 

so that X^ n '{t) converges to X(t) in probability. 

Remark 3.5. If Err(X(t), X^>(t)) tends to as n goes to infinity for all t G [-T,T} C - 
then 



{^ (n) (*)}te[-T,T]9 {X(t)} t< z[-T,T] 



f-d. 



f-d 



where —* denotes convergence in distributions of all finite dimensional marginals. 

Proof. We first state a lemma which is an implication of the inequality 

(x + y) p < x p + y p , 0<p<l, x,y>0. 

Lemma 3.6. Let < p < 1 and G L p (R d ) for i = l,...,n with n G N. Then 



i=l 



< En/*- 

LP i=l 



Now fix any ti, t m G [— T, T] 9 and Ai, X m G M. If < a < 1, we get 

' / m m 

Err a [Y,\X(t 3 )^^X {n \h, 
\i=i i=i 



(n) 



3=1 

m 



3=1 



m 

< E M 



(n) 



j'=i 



n — > oo. 



For 1 < a < 2, we can use Minkowski's inequality and get 

/ m m 



3=1 



3=1 
m 



3=1 



m 

<Ew 

L« 3=1 



f{n) 



3=1 



n — > oo. 
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Analogously for a = 1, 

(m m \ 

3=1 i=i / 

Therefore, for any Ai, A m G M we have 

m m 
3=1 3=1 

in probability which implies the convergence of all finite-dimensional marginal distri- 
butions. □ 



3.1.2. Poisson random measures. 

Let $ be a Poisson random measure with intensity measure and \& be the random 
sequence of the Poisson point process corresponding to the Poisson random measure. 
Furthermore, assume that f t is measurable on [— A, A] d for each t G [—T,T] q . Then, 
by the Campbell theorem (see [T8] , p. 103), we have 



£rri(X(*),! (n) (t)) 



f t (x)-fT\x) Q(dx) 



-A,A] d 
( 



E 



f t {x) - ~f\ n \x) 



\[~A,A] d 



*(*)-/ t (n) (x) 



> E 



E 
E 



/ t (x)$(dz) - / / t w (z)^W 

-A,A] d J[-A,A] d 

X{t)-X {n) {t) 



that is we can control the mean error between X(t) and X^it) in the L 1 -sense by 
finding error bounds for Ern{X(t),X^{t)). 



3.1.3. Exploiting the spot variable representation of the second moment of X(t). 
We assume that the second moment of the random field exists and recall formula © , 
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p. El 



E (X(t) 2 ) = / fi{y)VM{L>{y))\{dy) + / f t (y)E(L'(y))X(dy) 

J[-A,A] d \J[~A,A] d 



which implies 



E{X{t) -X (n) (t)) 2 



< 



(ft(y) -fT\y)) 2 v™(L'( y ))\(dy) 

l-A,A] d 

( r V 

+ J {m-R n \y)ML\y))Kdy) 

\[-A,A] d J 

(ft(y)-fT\y)) 2 V™(L'(y))\(dy) 

-A,A] d 

(My) ~ ft\y)?Kdy) ■ J (E(z%))) 2 X(d y ) 

-A,A] d l-A,A} d 



where we used the Cauchy-Schwarz inequality in the last inequality. If 



Var(L'(y)) < ci < oo, \/y E R d and f (E(L'(y))f \{dy) := c 2 < oo, 

J[-A,A] d 

then we get 

E(X(t) -X^\t)f) 112 < ( Cl + C2 )V2|| /t _ gt \\ L2 = ( Cl + C2 )^Err 2 (X(t),X( n \t)). 



Example 3.7. We choose again the Lebesgue measure for v and the characteristic 
triplet (a, 0,U) from Example 12.51 p. M We have 



Var(L'(y)) 



2 



-: Ci and 



(E(L>(y))f\(dy) 



(2AY 

e 2 



c 2 , 



[—A,A] d 



such that 



(E(I(t)-lW(t)f) 1/2 < i (1 + (2A) d ) 1/2 Err 2 (X(t),X (n \t)). 
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3.2. Step function approximation. 

For any natural number n > 1 and k = (ki, kj) £ 7L d with — n < k\, kj < n, let 



A fc 



.4 



n 



A 



n 



fci — , (&4 + 1)— ] x • • ■ x 

n n 



fed—, (fed + 1) — 
n n 



We define the step function 



|fc|<n 



to approximate / t , where |fc| < n is meant to be componentwise, i. e. — n < ki < n 
for i — 1, d Then we have 

X^\t)= J f t {n \x)A(dx)=J2ft^)A(A k ). (8) 

[_ A ^ ]d \k\<n 



The following theorem provides error bounds for Err a (X(t), X {n) (t)) for Holder- 
continuous functions f t . 



Theorem 3.8. Assume that < s < 2, t/ie control measure X is the Lebesgue measure 
and the functions f t are Holder-continuous for all t £ [— T, T] 9 , i e. 

|/t(^)-/ t (2/)|<a-|N-2/||?, MM", t£[-T,T]« 

for some < 7 t < 1 anc? C 4 > 0, where \\ • \\2 denotes the Euclidean norm. Then for 
any t £ [— T, T] 9 we have for all n > 1 i/iai 



(9) 



Proof. Since 



-IW(t) = / f/ t (x) - m)ZA k (x)) A(dx) 
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we have 



(Err a (X{t),X<- n \t)) 



-A,A] d 



-A,A] d 



/t(*) " £ /t(£*)W*) 
\k\<n 

|fc|<n 



dx 



x) 



dx. 



For each x G [— A, A] d , there exists exactly one k = k(x) with \k\ < n and x G A^. 
Hence 



|fc|<n 



K/t(x)-/ t (e fc ))i s i Ai (x) 

£ !/*(*)- /*(&)!* *A fc < 

|fc|<n 



.1- ) 



which implies 

'Err s (X(t),xW(t)) 



E (/*(*) - /*(&)) *a» 

|fe|<n 

= E / i(/*w-/*(e fc ))r<fc 
< a e/ ii^-aiir^ 



fix 



l fc l<™A, 



A/n A/n 



a£ /•••/ M+'-'+^^dw-dyi. 



I fc l<™ 



As ( 7t s)/2 < 1, we have (j/? + ■ • • + ^) {7t ' )/2 < y? s + ■■■ + yj ts and hence 

A/n A/n 

J - J (y 2 i + --- + ylf ts)/2 dy d ---d yi 







11 



A/n 



- I / y? s dyi 



av^ 1 r «... d MV +7ts 







7 t s + 1 V«/ 



(10) 
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Therefore, we get 

Err s (X(t),X^(t)) < 
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2 d C t d 



1 + Jts 
□ 



l/s 



i \ ft 
J^t+d/s I i_ 



n 



□ 



Remark 3.9. It suffices to consider a control measure A proportional to the Lebesgue 
measure (cf. Example 12.51 p. 12.51) . that is 

X(drj) = c ■ drj, c > 0. 

In this case, one has to multiply the upper bound in by c 1 / 2 . If the control measure 
A is not the Lebesgue measure, then, in general, the integral 

A/n A/n 

(y 2 1 + --- + y^ ts)/2 X(d(y 1 ,..., yd )) 



in (TTOj) cannot be calculated explicitely such that one would have to include it in the 
upper bound of the approximation error. 

Remark 3.10. In the proof, one can estimate 

A/n A/n 



/••■/ (y! + --- + yl) M/2 dy d ---d yi 



alternatively by 



A/n A/n 



y\\2 dy d ---dyt < 



1 



\\vh<$Vd 

and calculate the last integral by using polar coordinates. This yields 



Err s (X(t),X^(t))< 



where 



D(d,s) := < 



' 2 d C t ^ tS+d ^^ 
Its + d 



2 i/ S) 

d_ 3 r(3/2) 



l/s 



jpt+d/s 



11 



It 



D(d, s), 



7T 



r(d/2) 



l/s 



^d-7/2 T(3/2) 

71 ' r((d-i)/2) 



l/s 



d = 2, 
d = 3, 
d = 4, 

d > 5 odd, 
d > 5 even. 
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It is straightforward to show that for d = 2, this estimate is worse than the one in 
Theorem 13.81 if 

2 v ' ! 7t s + d ~ 

This is illustrated in Figure CD For d > 3, however, the estimate in Theorem 13.81 

always performs better. 




Figure 1. Combinations of (■y t , s) for d = 2. For all combinations in 
the dark grey area, the constant in Theorem 13.81 performs better, for 
all other combinations (in the light grey area) it is vice versa. 



Remark 3.11. Suppose that the conditions of Theorem 13.81 hold true. If the support 
of ft is not compact, we first need to estimate 

X(t) = [ f t (x)A(dx) 



by 

X K (t) = J f t (x)A(dx) 

[—K,K] d 

For K > large enough, the approximation error is small since 



Err s {X{t),X K {t)) 



\ft{x)\ s dx 



-K,K}° 



0, K 



/ 
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and f t Ji n) G L s ([-A,A] d ,\). Let e > 0. If 1 < s < 2, choose K > such that 
Err s (X(t),X K (t)) < e/2. We can apply Theorem EH to X^(-) such that 

Err s (X K (t),XP(t))< £ - 

for n G N large enough. Then 

Err s (X(t),X%\t)) < Err s (X(t),X K (t)) + Err s {X K (t),X%\t)) < £ + £ = e . 

If < s < 1, choose K > such that Err s (X(t), Xjc(t)) < e s /2. Again, we can apply 
Theorem 13.81 to Xk(-) such that 

Err s (X K (t),xP(t))< £ ^ 

for n G N large enough. Then 

Err a {X{t),X$\t)) < (Err a (X(t),X K (t)) + Err s {X K {t),X%\t))) llS 




Remark 3.12. Theorem 13.81 provides a pointwise estimate of the approximation error 
for each t G [—T,T] g . We can obtain a uniform error bound as follows. 

Assume that 7 := inf 74 > 0. Then for each t G [— T, T] q , f t is Holder- 

te[~T,T]i 

continuous with parameters 7 and some constant C 4 * > 0. Set C := sup C t *. 

te[,T,T] d 

Then Err s (X(t), X^^t)) can be estimated by ([9]) with C t and 7 t replaced by C and 
7- 

One can also consider the integrated error 

Err s (X,X {n) ) := y Err 3 (X(t), X {n \t))dt 

[-T,T}i 

and multiply the error bound by (2T) q . 

Remark 3.13. Assume that < s < 2 and the functions f t are differentiable with 
||V/ t (x)|| 2 < C t for all x G [-A, A} d and t G [— T, T] 9 . Then for any t G [— T, T] 9 , © 
holds for all n > 1 with 74 = 1, that is 

Err a (X(t), XW(t)) < ( • - 

\ 1 + s ) n 

since f t is Holder-continuous with C f and ^ t = 1. 
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3.3. Approximation by wavelet series. 

3.3.1. Series representation of kernel functions. 

Let s>0,/ t e L s (R d ), t G ffi 9 , and let {&} i& r be a basis for L s (R d ), where I is an 
index set. Then f t can be represented as 

/t = J>i'6 ( n ) 

for certain constants a* G K. In order to approximate f t , one can truncate ( TTTT ) 
such that it consists only of a finite number of summands. In [I], the trigonometric 
system is used to approximate the kernel function of certain stable random fields. In 
this paper, we will go another way and analyse whether a wavelet system may be 
appropriate for the simulation of random fields with an infinitely divisible random 
measure as integrator. 

3.3.2. Haar wavelets. 

In this paper, we use the so-called Haar basis to approximate the kernel functions. 
For a detailed introduction into wavelets, see for example [20], [5] and [4]. 

Definition 3.14. Consider the function 

and the corresponding mother wavelet defined by 

tfHaar^ . = ^Haar Q x + A) - <^ Haar (2x - A), X G R. 

The resulting basis H := {</? Haar } U {^J^l k } k^N , is called Haar basis for 

2 k <j<2 k+1 -l 

L 2 ([-A, A}), where #gf r (x) := 2 fc / 2 * Haar (2 fc (x + A) - (1 + 2 • 

Let \1/ := y9 Haar , vl/ 1 := v|/ Haar and E be the set of nonzero vertices of the unit cube 
[0, l] d . Consider the multivariate functions \l/ e , e = (ei, ...,ed) G E, defined by 

* e (xi, x d ) := * ei (^i) • • • ^"(xd), x G M d . 

Let x = (xi,...,Xd) and a = (A, A) T , c = (1, 1) T G M d . Translation by j = 
(./: 3d) and dilation by 2 k yields * e „ 2fcc>fc (x) := 2 kd ^ e (2 k (x - a) - A{c + 2j)), 
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2 fc < ji < 2 fc+1 — 1, k G N , z = 1, ...,d, e 6 £, that form an orthonormal basis of 
L 2 ([-A, A] d ). Then, each / G L 2 ([-A, has the expansion 



/* = (/*,**)**+ EE E (/t^-^*)^-^. (12) 

eG-B fc=0 2 fc <ii<2 fc+1 -l 
i=l,- ,d 

in the sense of £ 2 ([— A, A] d ) convergence, where 



:= (2l)^' ^H 1 ^ ( 13 ) 
(/*.*J-*c*) : = / f(*)*j-»c*(*)ta- (I 4 ) 

It can be shown that any function / G L max { s '?>}([— A, A] d ) for some p > 1 can be 
represented by a wavelet series of the form ( IT2l) in the sense of L max ^' p >([-A, A] d ) 
convergence. However, there exist examples of functions for which (fT2l) does not hold 
in particular for s = p = 1, cf. [5], p. 7. Therefore, we restrict our setting to kernel 
functions / t G L max ^ p >([-A, p > 1. 



As noted, we want to use the expansion (JT2J) in order to approximate the kernel func- 
tions ft by truncating the (potentially) infinite sum to a finite number of summands. 
The goal is then to find an upper bound for the approximation error. 

3.3.3. Approximation by cutting off at a certain detail level. 

Consider a kernel function f t G L max ^ s ' p ^([— A, A] d ) for some p > 1 with corresponding 
Haar series 

oo 

/* = (/*, EE E (/*» *u*c, k )yu*c, k - 

e£E k=0 2 k <ji<2 k+1 -l 
i=l,- ,d 

The idea is now to cut off this series at a certain detail level k = n, that is to 
approximate the kernel function f t by 

n 

%L = (/*, + E E E (/* *U*c, k )*u*c, k - 

e&E k=0 2 k <ji<2 k+1 ~l 
i=l,'" ,d 

The following lemma provides an upper bound for the approximation error of bounded 
kernel functions by applying the cut-off truncation method. 
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Lemma 3.15. Let s > and d > s. Assume that M t := sup < oo. Then 

xe[-A,A] d 

for n6N 



(fe) VS • ■ M t ■ (2A) d /° ■ {^) n , < s < 1, 
^- i .d-M t -(2A) d / s -(^ T ) n , 8 >1. 



Wft-L 



(n) I, 
t,cut\\L s 



< 



Proof. Let s > 1. We have 



ll/t-/ 



(») II 



eS-B k=n+l 2 fc <j i <2 fc+1 -l 
i=l,- ,d 



j — 2 k c,k 



< 



E E E 

egE fc=n+l 2 fe <j i <2 fc+1 -l 
i=l,- ,d 



j-2 k c,kJ\ \\^ j-2 k c,k\\L a 



Now |(/t,^_ 2 fc c ,fc)l can be estimated by 

l(/*.*J-*c,*)l < M * / l«JW af )l daf = M * 



2kd/2 

(2A) d / 2 



(15) 



■ 2- fed • (2A) f 



and 



^J— 2 k c,k 



is equal to 



I j-2 fe c,fc|lj> 



|^ e _ 2 , c , fe (x)|Mx 



9 M 

Z 2 fed ... i 

2~"r • (2i4)- 



{2A) d /i 



Thus 



n/t - /?,iiu* < E E E M * • 2 " M/2 • 2¥ • 2 " ¥ • ( 2 ^) f 



eS-E fc=n+l 2 fe <ji<2 fe+1 -l 
i=l,- ,d 



(2 d -l)M t (2A) d/s J2 d-2 k -2~ kd/s 

k=n+l 

oc 

(2 d - 1) M t d{2A) d ' s (2 1_d/s ) fc • 



fc=n+l 
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Since d > s, we have 1 — d/s < and by using the geometric series formula, we get 



Now let < s < 1. By Lemma [3761 we have 



\\ft ftjutllh 



EE E (^-2^)^-2^ 

eS-B fc=n+l 2 k <ji<2 k+1 ~l 
i=l,- ,d 



L* 



< 



EE E 

ee_B fe=n+l 2 fe <ji<2 fc+1 -l 
t=l,- ,d 



I s V[/ e 

j-2 fc c,fc/l || j-2 fe c,fc|lxs 



By using the estimates for the wavelet coefficients and the L s -norms of the wavelets 
from above, we get 



\\ft-f&\\h < (2 d -i)d(M t (2A)^y ( 2 ^) fc 



k=n+l 



and finally 



(*-l) 



(2 d ~ s - 1) 



2 d/ 



s-l 



□ 



□ 



If we make further assumptions about the kernel function f t , we can improve the rate 
of convergence of the upper bound. 

Corollary 3.16. Assume that f t is Holder- continuous with parameters C t and ^ t for 
all t G [-T, T] g . Then for neN 



wit-filth* < < 



2 d -l N ^ 



T ) 3 ■d^/^-C t .(2A) d /^.(^)\ 0<s<l, 



2 d -l . jl+7t/2 . n ■ (9 A\d/s+lt . I 1 ^ 



S > 1. 



Proof. We estimate \ (f t , ^_ 2 fc cfc )| as in the preceeding lemma. Let 

B := {xe[-A,A] d :^_ 2 ^ k (x)>0}, 
C := {xe[-A,A} d :^_ 2kCtk (x)<0}. 
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Then we have 



lift, *l-2*c,fc) 



l-A,A] d 
2kd/2 



ft{xm_ 2kc jx)dx 



2kd/2 r 2 kd / 2 

B Mx) j^dx - J c f t (x) jMj&dz 



ft(x)dx - / f t (x)dx 

B JC 



2kd/2 

< - — ~ ttt max 



(2A) rf / 2 



ft(x)dx - max(f t (x))\C\ 



ft{x)dx - imn(f t (x))\C\ 



B 



x&C 



Here, \C\ denotes the volume of C. We now estimate the two quantities in the 
maximum. We have 



ft{x)dx - max(/ t (x))|C| 



' B 

< max 



max(f t (x))\B\ - max(/ t (x))|C| 

x£B xSG 



m Mft(x))\B\ 

x£B 



max(Mx))\C\ 

x£C 



and 



ft(x)dx - min(/ t (x))|C| 



' B 

< max 



xec 



m a x(f t (x))\B\ - mm(f t (x))\C\ 

xeB xeC 



mm(f t (x))\B\ 

x£B 



mm(f t (x))\C\ 

xec 



Therefore we get 



< 



I 

2fcd/2 

(2A) d l 2 



2 k c,k) 



max 



max(f t (x))\B\ - max(f t (x))\C\ 

x&B x&C 



max(f t (x))\B\-min(f t (x))\C\ 

xGB xdC 



min(/,(x))|S|-max(/ t (x))|q 

x&B x£C 



mm(f t (x))\B\-mm(f t (x))\C\ 

xeB xGC 



Let X\ G B and x 2 G C such that the maximum in the last inequality is attained. 
Furthermore, we have 



\B\ = \C\ = --(2A) d .2 



-kd 
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Then, since f t is Holder-continuous, we get 



j-2 k c,k) 



9W/2 n kd / 2 1 

< j2Ayf2 \M*i) ■ \*\ - f*M ■ \C\\ = j^j- 2 ■ -i2A) d ■ 2- M |/*(x 1 ) - f t (x 2 )\ 

nkd/2 1 nkd/2 1 , . ™ 

_ ^d/2+7t(7 4 2 d / 2 +7t-l c /7t/2 



2^/2+7* 

The remainder of the proof is analogous to the one of Lemma [3.15L □ □ 

Corollary 3.17. Assume that f t is differentiable with ||V/t(x)||2 < C t for all 
x e [-A, A} d , C t >0 andte [-T, T} q . Then for n e N 

1/8 



< Jlfafe) ■d^.C t .(2Ay/^.(^_) n , 0<s<l 
U ~ kcUth '-<J^/**.C t .(2A)«*.( M ", s>l. 



3.3.4. Near best n-term approximation. 

Taking a wavelet basis for L max ^ s,p ^([— A, y4] d ), p > 1, has advantages in particular in 
the representation of functions with discontinuities and sharp peaks, that is functions 
with a certain local behavior. By simply cutting of at a certain detail level, this 
advantage is not honored. In view of (fl5l) . we may expect that we have to calculate 
less Haar coefficients if we approximate the kernel function f t by a truncated Haar 
series that contains those n summands (/, ^ e _ 2 k c k)^'j-2 k c k w ^h the largest val- 
ues ||(/, ^_ 2 fc c ,J^_2fc c ,JU s = l(/'^-2*c,fc)l ^ e j~2><c,k ■ An approach to use such 
a truncation in order to approximate functions is presented in [4j, p. 114 ff. We 
summarize the main statements. 



Consider a function S defined by 

S= a l^U^ S e ,* eM > V(e,j,fc)eH (17) 

(e,fc,i)GH 

where E := {(e, k,j):eeE,ke N , 2 fc < ji < 2 k+1 - l,i = 1, d) and #H < n for 
some neN. Hence, (fTTjl is a linear combination of n Haar wavelets. We denote by 
E„ the set of all the functions S defined as in (1171) and let 



a s n (ft) := inf ll/-^^ 
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Now, we truncate the wavelet expansion ([T2 



for which the absolute value of |(/, \P 



j-2 k c,k> 



of ft by taking those n summands 
is largest and denote the 



j-2 k c,k 



truncated sum by fl . In [19j, the following theorem was proven which shows that 
this truncation is a near best n-term approximation. 

Theorem 3.18. Let 1 < s < oo. Then for any f G L s {[-A,A] d ) we have 

\\f-f(-)\\ Ls < Cl (s,d,A)a s n (f t ) 
with a constant C\(p,d,A) > only depending on s, d and A. 



If the sequence ||(/, ^_ 2 ' 
< r < oo, that is 



c,k' 



\l> e u 

j—2 k c,k 



> is in the Lorentz space wl 7 

1 J j,k 



#{ 0;fc,e):{||(/,^ fe )^J is >£}<(- I , V £ >0, M>0, (18) 

with a certain additional condition on r, then cr^(ft) can be bounded from above as 
shown in [4], p. 116. 

Theorem 3.19. Let 1 < s < oo and f G L s ([-A, A] d ). Furthermore, let 

and u > wit/z 1/r = u + 1/s. Then 

a s n (f)<C 2 (s,d,A)M(^j , n = l,2,-.-, 

w^/i a constant C^s, d, A) > on/y depending on s, d and A and M being a constant 
satisfying (jTM . 



Combining Theorem 13.181 and Theorem 13. 191 yields an upper bound of the approxima- 
tion error by using the near best n-term approximation with a rate of convergence of 
O ((l/n) s ). In [TO], we obtained the following formulas for C\(s, d, A) and 6*2(5, rf, A) 
in the case that / G L s ([0, l] '): 



d(M) 
<7 2 (M) 



2+1-2" 



-2 



(2 d - 1) 



max I s, 



s - 1 



2 (2 T / S - 1) 
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We note that on the one hand, these constants are not sharp and may be quite large, 
and on the other hand, we would have to find a value of M in (fTHl) as small as possible 
for each kernel function ft or for certain classes of kernel functions, which is not so 
easy to determine. Therefore, we suggest an approach which is not based on the error 
estimate with those constants, but still determines an approximation at least close to 
the near best n-term approximation while keeping the desired level of accuracy. 

It is clear that || (/, i &j_ 2 k c k)^'j-2 k c JU S § oes *° zero as the detail level k goes to 
infinity. This means that the n largest values || (/, s &j_ 2 k c k)^j-2 k c k\\ L " are likely to 
be found for small values of k. 

We now assume that s > 1, d > s and M t := sup \ft{%)\ < °° and choose e > 

x£[-A,A] d 

as the desired level of accuracy. The following derivations are analogous for the case 
< s < 1. 

From Lemma 13.151 and its proof, we get 

oo 

\\ft-ll%t\w < EE E \\u^u^u*c, k \w 

e£E k=n+l 2 fe <j 4 <2 fe+1 -l 
i=l,- ,d 

if and only if 



\n(e(2 d / s ~ 1 - 1)) - ln((2 d - l)dM t {2A) d l s ) 
U ~ - ln(2 d / s - 1 ) 

We take 



m t :-- 



\n(e(2 d / s ~ 1 - 1)) - ln((2 d - l)dM t {2A) d l s ) 



_l n (2^-i) 

where \x] is the integral part of x, as our minimal detail level to obtain the desired 
level of accuracy and add / G No detail levels to the truncated wavelet series at detail 
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level m t , that is we consider 



Amt+l) 
J t,cut 



mt+l 



(/*,**)** + X)E E u^u^u^h- (19) 



Furthermore, we define 



eS-B fc=0 2 fe <j i <2 fe+1 -l 
i=l,— ,ci 



mt+l 



c : = ii (/ t> + EE E ii(/*»«;-^*)«5 



(20) 



eGS fc=0 2 fc <ji<2 fc+1 -l 
i=l, ••■ ,<i 



and 



d ■.= \\(f t ,^m\ L s + J2 E E ik/^^^ju* 



eS-B k=mt+l 2 k <ji<2 k+1 -l 
t=l,.» ,d 



By Lemma [3. 151 the corresponding level of accuracy of (JT9|) is at least 



2 d - 1 

2 d/s-l _ I 



d-Mf (2A) 



d/s _ 



1 



2<i/s-l 



mt+l 



Now we take the n largest summands from (J20| . that is from 

{iiu^^^iU4U{ii(/*'^ c , fe )^ Clfe ii^k fe , 

and denote them by a\,...,a n . The remaining summands are denoted by a n+ \, a n+ 2, ■ 
and the corresponding summands from (fl9l) by bi,...,b n , b n+ i, b n+2 , ... . The number 
n G N is chosen to be the smallest number such that 

n 

C — E a i ^ e ~ e *f 



i=l 



We define 



fV :=E^ 



i=l 



which is close to the near best n-term approximation if I is chosen large enough since 

II [ftj ^ j-2 m t+ l ,m t +d^ j-2 m t+ l ,m t +l W LS 

goes to zero as I goes to infinity. 
Then we have 



II/-/ 



(n)|| 
t \\L S 



i=n+l 



< E CLi = C — E CLi + D < 6 — E* + E* = E. 
j^a i=n+l i=l 
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3.3.5. Implementation. 

When implementing the wavelet approach, one more problem has to be considered 
which we discuss now. 

Let / be the set of the indices (e,j,k) for which the summands (f t , ^|_ 2 fc c fc)^j_2 fc c k 
are part of the approximation . Then we can write 

if (f t , vj/*)i|)* is included in the truncated series or 

ft ) = E ^ e j-2 k c,k)^ e j-2 k c,k 

(e,k,j)el 

if it is not included. 

In order to approximate the random field X, we use 
or 

(e,k,j)ei { _ AA]d 
if, again, (f t , is not included in the truncated series. 

Since the Haar wavelets ^j_ 2 k c & are simple step functions, the integrals 

Vj-2»c,kM dX ) 



[—A,A] d 

can be easily simulated although they are not independent: Let z G N be the finest 
detail level of the wavelet approximation. Then all of these integrals can be built up 
from 

J A(dx) = A^ -A + k^,-A+{k + l)£) j , £; = 0,...,2 2+1 -l 

which are independent because the sets [—A + k-^, —A + (k + 1)^?) are disjoint. 

However, the wavelet coefficients (/ t , ^j_ 2 k ck ) cause problems if no closed formula of 
the integral of the kernel functions f t over cubes is known. In this case, they have to 
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be determined numerically by using the fast wavelet transform (see for instance |20j, 
pp. 134). 

This results in a further approximation error which we need to estimate. When the 
detail level at which the wavelet series is cut off is equal to n, the input vector of the 
fast wavelet transform consists of integrals of the form 

2(n+l)d/2 

where cube is a cube of side length 2~ d ( n+1 \ 

We now assume that we have calculated these integrals with a precision of 5 > and 
denote the wavelet coefficients computed by the fast wavelet transform by (f t ,ty*) 
and C/t,^_ 2 >= C)fe )- 

When applying the fast wavelet transform, the value of each integral is used 2 d — 1 
times at each detail level k to calculate the wavelet coefficients (ft,^ e _ 2 k ck ), e £ E. 
There are 2^ n+1 ^ d such integrals. 



When s > 1, the precision of 

n 

- j-2 k c,k> j-2 k c,k 



f& = (/*,*•)*• +EE E (/*> *u^i 



e£E k=0 2 fe <i i <2 fc + 1 -l 
£=!,••■ 4 



to approximate 



/S* = EE E u^u^u^k (2i) 



e&E k=0 2 k <ji<2 k+1 -l 
i=l,- 4 



IS 



EE E i(/*^c,J - (/t^i-^c,Jiii^ c>fc ii^ 



egE fe=0 2 fc <ji<2 fc+1 -l 
i=l,- 4 



(n A\d/s-d/2 n nkd/2 
n(n+i)d \ffV A _l_ (0^ 1"\ n(n+l)d r Icy a \d/s-d/2 c)kd/2-kd/s 

- z 2( n+1 ) d / 2 1 2(" +1 W 2 ^ ^ 



fc=0 



< f ( 2A)^^/ 2 2("+^/ 2 (i + (2* - i) 2 ^;::;'- 1 ) *, *>i, (221 

\(2A) d / 2 2( n+1 ) d / 2 (l + (2 d -l)(n + l))5, s = l. 
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When < s < 1, we can use Lemma [3761 and get 

\\Mt-ftlth* < (\(7^) - (ft,**)\ s \mi* 
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2 k c,k)\ S \\^ e j-2 k c,k\\L* 



e&E fc=0 2 k <ji<2 k+1 ~l 
i=l,- A 



/ 9 (ds~d)(n+l) _ i \ V* 

< (2A)d / S -d/2 2 (n + l)d/2 h + (2* - 1) 2d ,_ tf _ 1 ) & 

Let s > 0. We choose e\ > and e 2 > such that s x + £ 2 = e if s > 1 and ef + £3 = £ s 
if < s < 1. Furthermore, we choose the detail level n so large (by using the formulas 



in Section 13.3.31) such that 

We approximate the elements of the input vector for the fast wavelet transform with 
a precision of 



( 2A )d/ s -d/2 2 («+l)d/2 ^ 1 + ( 2 d_i) g^j^^llLll^ 



T77, 0<3<1, 



5 = < 



£2 



(2A) d / 2 2("+ 1 ) d / 2 (l+(2 d -l)(n+l)) ' 

£2 

(2A)"/'-''/22("+i)''/^l+(2''-l)2 ( ' l ~^" + ^- 1 J' 

Then we have for s > 1 



s = l, 

s > 1. 



Il/t ~~ /i.eiilU 3 — ||/t — /i.cut II L° + Wft Jut ~ ft,cLt\\L s 1 

and for < s < 1 



?(n) 



?(»*) ?w 



ll/t — AeutlU 8 — ( ll/t — /i,ciitll£ s + ll/t,cut — ft,cutWl s ) 



?(n) 



(n) ?(n) 



We summarize this result in the following algorithm. 



Algorithm 

Let Mt := sup |/t(^)| < 00 and d > s. Choose e > as the desired level of 

xe[-A,A] d 

accuracy. Choose £i,£ 2 > such that e = e± + e 2 if s > 1 and e = (e\ + e^) 1 ^ if 
< s < 1. 
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(a) Let 



m t :-- 



ln{e 1 (2 d - s -l) 1 / a )-ln{(2 d -l) 1 / a d 1 / a PIt(2A)-) 
(1-d/s) ln(2) 

ln(ei(2 d / s - 1 -l))-ln((2 d -l)dA^(2A) d / s ) 
(1-d/s) ln(2) 



< S < 1, 
S > I. 



and choose a number I G N that increases the detail level m t . 
(b) Calculate the wavelet coefficients for 

mt+l 

(f u + E E E ^u^u^ 

e£E k=0 2 k <ji<2 k+1 -l 
t=l,— ,d 



using the fast wavelet transform with a precision of 



5 = < 



£2 



( 2j4 )d/ s -d/2 2 (n+l)d/2 ^ 1 + ( 2 d_l) gjj^g^llLll^ 



£2 



(2j4) d2(n+l)d/2[ l+( 2 d-l) 2 t-d/at. 1 ! 



: C2 

( 2A )d/ s -d/2 2 (n+l)d/2^ 1 + ( 2 d^ 1 )2^^^+|)-iy 



< s < 1, 
s = l, 

s > 1. 



(c) Take the n largest summands from 

m t +l 

C := \\(f^)**\\ L s + E E E IMiV^jVj-vJl" 

e£E k=0 2 k <ji<2 k+1 -l 
t=l,- ,d 

and denote them by ai,...,a n . The corresponding summands from 

mt+l 

(7^*w + E E E (/*.^c*)«i-^ 

ee£ fc=0 2 fe <i i <2 t; + 1 -l 
i=l,- ,d 

are denoted by bi,...,b n . Choose the number n such that 

n 

c - E°* - £i ~ £ *' 

where 

' (.l^) 175 " " M t ■ (2A) d/s ■ (^r +l , < 8 < 1, 
^L T .rf-M,.(2A)^.(^ T ) mt+Z , S >1. 



(d) Take / t = X/ILi 6j as the approximation for f t . 
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Remark 3.20. (a) Assume that d EN and f t is Holder-continuous with parameters 
C t and 7t for t E [— T, T] q . Then the algorithm can be applied with m t and e\ 
replaced by 



mt 



ln(2£i(2 d + a Tt-l) 1 / s )^ln((2 d ^l) 1 / s d 1 / a +Tt/( 2a )Ct(2A) d / a +Tt) " 
-ln(2 d / 3 +T«) 

ln(ei(2 d / a +Tt+ 1 -2))-ln((2 d -l)d 1 +Tt/ 2 C t (2A) d / a +Tt) ' 
-ln(2 d / 3 +Tt) 

1/S -d^.C t -(2A)^.(^r +l 

2 d/s+- lt +i_ 2 u °i l2 d / s +^/ ' 




2 a -l 
J V 2 d+s T«-l 

2 d -l 



< s < 1, 

s > 1, 

< s < 1, 
s > 1. 



(b) Assume that d G N and / t is differentiable with ||V/ t (a:)||2 < C t for all x E 
supp(/t) with C t > and t E [— T, T] 9 . Then the algorithm can be applied 
with m t and e* replaced by 

!+s iM/s\ i_/^nd 1 \l/s j3/f2s1^> /r, /l\d/s + l\ ~1 

< s < 1, 



ln(2ei(2 d + 3 -l) 1 / s )-ln((2 d -l) 1 / 3 d 3 /( 23 )C i (2A) d / a + 1 ) 

-ln(2 d / s + 1 ) 
ln(ei(2 d / a + 2 -2))-ln((2 d -l)d 3 / 2 C t (2A) d / 3 + 1 ) ' 
-ln(2 d / s +l) 

1/ S 



s > 1, 

mt+l 



2 
2 

2 d / 



mt+l 



, < s < 1, 

s > 1. 



We conclude this section with the main result. 

Theorem 3.21. Assume that s>0andf t E L ma ^ s ^([-A, A} d ) for some p > 1. Lei 
X be an infinitely divisible random field and the control measure X of the infinitely 
divisible random measure be the Lebesgue measure. Let e > 0. If ft is calculated 
using the algorithm mentioned above, then 

Err s {X{t),X in \t)) < e, Vt G [— T, T] q . 
4. Simulation study 



For the simulation study, we used two different types of kernel functions for a-stable 
random fields of dimension d = 2. The first one is an Epanechnikov-type kernel 
function defined by 



ft(x) 




t 



||x - t|| 2 < a, 
otherwise, 



(23) 
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where a > and b > 0, whereas for the second one, we take 

f( tl ,t 2 )(xi,x 2 ) = b(a - \x\ - h\){a - \x 2 - t 2 \) 

■ l{a-|xi-ti|>0, a-\x2-t 2 \>0}( x li x 2) 

where a > and b > 0. Examples of both kernel functions are plotted in Figure [2l 



(24) 





Figure 2. The Epanechnikov-type kernel function ( |23l) (left) and the 
kernel function ( 1241) (right). 

The main difference between these two types of kernel functions is that one can derive 
a simple formula for the integral of (1241) over squares, but not for the integral of (T23l) . 
This does not affect the step function approach since the kernel functions are only 
evaluated there at the points but it does affect the wavelet approach because the 
input vector consists of such integrals of (123]) and (T24l) over squares. Therefore, we 
have to expect a loss in computational performance for kernel (T23l) with the wavelet 
approach in this case. 

Both functions (1231) and ( 1241) are Holder-continuous with parameters (Ci,7i) = 
(2ab, 1) and (62,72) = (V2ab, 1), respectively. We fixed a = 1.5, (3 = and 
[-T, T] 2 = [—1, l] 2 for both types of kernel functions. 

For the remaining parameters, we started with the following configuration: b = 1, 
a = 1 and e — 1. Furthermore, we divided [—1, l] 2 into an equidistant grid of 50 x 50 
points and chose I = for the number of detail levels to be increased. Two realisations 
of the 1-stable random field X with kernels f|23l) and ff24l) are shown in Figure [3l 

First, we kept all parameters fixed and determined the computational time depending 
on the number of realisations. For the step function approach, each realisation needs 
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Figure 3. Two realisations of stable random fields with kernel (l23l 
(left) and kernel ([24]) (right). 



the same computational time. For the wavelet approach, however, the wavelet coeffi- 
cients only have to be calculated for the first realisation and can be stored afterwards. 
Therefore, any further realisation needs less computational time. Table CD shows the 
results for both the Epanechnikov-type kernel ([231) and kernel (1241 ). By trial-and- 
error, we figured out that a combination of e = S\ + Ei = 0.99 + 0.01 performs quite 
good for the corresponding parameters in the wavelet algorithm in this case. 

Table 1. Computational time (in msec) for the first and further realisations. 





Kernel ([231 


Kernel lj24|) 


Step function 


28.5 


18.6 


approach 






Wavelet approach 


5337.0 


1090.0 


(first realisation) 






Wavelet approach 


13.5 


13.2 


(further realisations) 







Second, we focused on kernel ([241) and the computational time of any further reali- 
sation except the first one and varied subsequently one of the parameters a, m (the 
number of pixels per row) and e while all the other parameters were kept fixed. It 
turned out that the computational time decreased for the wavelet approach and in- 
creased for the step function approach when a was decreased. The computational 
time was equal for both approaches at about a = 1.8. For decreasing e, the computa- 
tional time for the step function approach increased much faster than the one for the 
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wavelet approach. Varying m affected the computational time of both approaches in 
a similar manner. The above results imply that neither of the two approaches out- 
performs the other one. For some combinations of the parameters, the step function 
approach was faster than the wavelet appraoch, for others it was slower. 

Finally, we increased the parameter I successively for a field with 10 x 10 pixels while 
all other parameters were kept fixed and investigated the computational time for the 
wavelet approach for any further realisation except the first one. Table [2] shows the 
corresponding results. 

Table 2. Computational time (in msec) for different values of I (kernel ( |24l) ). 

I 12 3 4 

Computational time 25.5 45.5 246.4 1044.8 4212.0 

One might have expected that the computational time tends to decrease if I is in- 
creased since the wavelet series usually consists of less summands when keeping the 
same level of precision. At the same time, however, more stable random variable 
simulations have to be performed for the calculation of the integrals 

/ vf e m+/M (dx). 

J[~A,A] d 

That is why for larger values of I, the computational time increases sharply. 

For a = 2, the random measure M is a Gaussian random measure and the random 
field 

X{t) = [ f t (x)M(dx) (25) 
Jr 2 

becomes a Gaussian random field. Therefore, we can use the simulation methods to 
simulate such fields, too. 

The kernel function 

f t (x) = b-e — , teR 2 

corresponds to the isotropic and stationary covariance function 
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:= Cov(X(0,0),X(/i,0)) = 2 / 6e"^ ■ be s a ete = rrafeV^, (26) 

cf. [H], P- 128. 

We now want to investigate how our simulation methods perform compared to the 
circulant embedding method (see [21]) to simulate Gaussian random fields with the 
given covariance function (1261 ) for a = 0.05 and 6=1 and [— T, T] 2 = [—0.5, 0.5] 2 for a 
grid of 100 x 100 points. We choose the circulant embedding method since it is exact in 
principle, and if exact simulation takes too much computational time, approximation 
techniques exist such that at least the one-dimensional marginal distributions are 
exact in principle. 

For the comparison of the circulant embedding method with the step function ap- 
proach and the wavelet approach, we simulate 1000 fields and estimate their mean, 
their variance and their covariance function for the distances 0, 0.01, 0.02, 0.5. 
We then compare the values to the theoretical ones and require that the estimated 
values do not differ from the theoretical ones more than 0.01. If at least one value 
differs more than 0.01, the precision level is increased. The following table shows the 
computational time for each of the three methods. 



Table 3. Computational time (in msec) for the circulant embedding 
method and the step function and the wavelet approach. 



Circulant embedding 


Step function approach 


Wavelet approach 


Computational time 48.43 


9588.29 


505.28 



As one can see, the wavelet approach and the step function approach take much more 
computational time than the circulant embedding method. This comes from the 
extensive calculations in the numerical integration of the stochastic integral (l25l) . For 
Gaussian random fields, it is therefore advisable to use existing simulation methods 
such as the circulant embedding method that exploit the specific structure of these 
fields. 
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5. Summary 

We presented two approaches to simulate a-stable random fields that are based on 
approximating the kernel function by a step function and by a wavelet series. For both 
approaches, we derived estimates for the approximation error 

In the simulation study we saw that for the first realisation of an a-stable random 
field, the step function approach performs better than the wavelet approach due to 
the initial calculation of the wavelet coefficients. For any further realisation, however, 
the wavelet approach outperforms the step function approach. 

Let us compare the rates of convergence of the step function approach and the wavelet 
approach more generally. If ft is Holder-continuous, the error estimate for the step 
function approach is 



for a constant Ci(d, C t , jt, s, A) > 0. In the simulation study, we have seen that 
increasing the parameter / in order to get closer to the best n-term approximation 
is not so advantageous. Therefore, we consider the rate of convergence for the cut 
wavelet series with error estimate 



with a constant C 2 (<i, C t , 7 4 , s, A) > 0. We note that we cannot compare the error 
estimates directly because for the step function approach, n determines the number 
of cubes ((2n) d ) that form a partition of [— while for the wavelet approach, 
n is the detail level. Therefore, we express the error bounds in terms of the number 
of summands of the step function approximation ((8)) of the random field, p. [03, and 
the wavelet approximation (I2T1 ). p. [301 respectively. In Table 31 we see that the rates 
of convergence in terms of the number of summands do not distinguish substantially 
between the step function approach and the wavelet approach. 

Let us consider the number of random variables that need to be simulated for a single 
realisation of a random field. For the step function approach, we need to simulate 
(2n) d random variables. For the wavelet approach, the number of random variables to 
simulate is equal to the number of cubes that form a partition of [—A, A] d in the finest 




(27) 




(28) 
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TABLE 4. Number of summands of (jSJ) and (l2lj) and error bounds for 
(1271) and (l28l) in terms of the number of summands. 





Step function approach 


Wavelet approach 


Number of summands 


u = (2n) d 


u = 1 + 2 d - 1 d(2' l+1 - 1) 


Error bounds 




o(( i 



detail level n: 2 d ^ n+l \ In the examples of the simulation study, much less random 
variables had to be simulated for the wavelet approach than for the step function 
approach which was a reason for the good performance of the wavelet approach. 

We want to make two remarks about the wavelet approach. First, we have seen 
that one drawback is that the computation of the input vector for the fast wavelet 
transfrom may take quite a long time if no formula for the integrals 

f t (x)dx 

is known, where C is a cube in M. d . In general, for an arbitrary wavelet basis {\E'*}j e /, 
we would have to calculate 

/ f t (x)^i(x)dx 

for some i & I which, in many cases, also requires numerical integration. 

Interpolatory wavelet bases can remedy this disadvantage since for this kind of 
wavelets bases, the wavelet coefficients basically reduce to evaluating the kernel func- 
tion at a certain point. However, the interpolatory wavelets themselves are no step 
functions any more such that the simulation of the integrals 

J[-A,A] d 

where is an interpolatory wavelet function, is much more complicated than for the 
Haar basis. 

Second, one could use adaptive wavelet methods in order to calculate the wavelet 
coefficients. This might decrease the computational time for the first random field 
realisation. However, in the simulation study we have seen that increasing the pa- 
rameter / has little advantage over the cut wavelet series (/ = 0) since the negative 



c 
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effect of the increasing detail level and thus the need of more stable random variable 
simulations dominates the positive one of less summands in the wavelet series. 
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